#This script will determine the residual value of the relationship between ES and a species abundance (tau=.975)
#Clear workspace

rm(list=ls())

in.dir = "C:/R/In/"
out.dir = "C:/R/Out/"

#Set working directory

setwd(in.dir)

#load necessary libraries

library("SDMTools")
library("quantreg")
library("adehabitat")

#Import the dataset for analysis

samples= read.csv("samplesmaster.csv",header=T)
georef = read.csv("georefmaster.csv", header=T)
s.fields = names(samples)
g.fields = names(georef)

#Perform quantile regression for each of 4 species against it's own ES values

cr.rq=rq(georef$cr_avgabund~georef$cr_es,data=georef,tau=.95,method="br",model=T)
sb.rq=rq(georef$sb_avgabund~georef$sb_es,data=georef,tau=.975,method="br",model=T)
gq.rq=rq(georef$gq_avgabund~georef$gq_es,data=georef,tau=.975,method="br",model=T)
lc.rq=rq(georef$lc_avgabund~georef$lc_es,data=georef,tau=.975,method="br",model=T)
sb2.rq=rq(georef$sb_avgabund~georef$sb_es+georef$lc_es,data=georef,tau=.975,method="br",model=T)
gq2.rq=rq(georef$gq_avgabund~georef$gq_es+georef$MODIS_mean+georef$lc_es,data=georef,tau=.975,method="br",model=T)
#Make a data frame from the regression residuals and the raw georef data
master=georef
master=data.frame(georef, cr.rq$residuals, gq.rq$residuals,lc.rq$residuals, sb.rq$residuals)

#Write out this data frame as a .csv file

write.csv(x=master, file=paste(out.dir,"quantregresidsbygeoref.csv",sep"/"),row.names=F)

#Perform MLRM with multi variables against avgabund

lc1.lm = lm(master$lc_avgabund~master$lc_es, data=master)
lc2.lm = lm(master$lc_avgabund~master$lc_es+master$bc05micro, data=master)
lc3.lm = lm(master$lc_avgabund~master$lc_es+master$bc05micro+master$lsc_mean, data=master)
lc4.lm = lm(master$lc_avgabund~master$lc_es+master$bc05micro+master$cr_avgabund, data=master)

sb1.lm = lm(master$sb_avgabund~master$sb_es, data=master)
sb2.lm = lm(master$sb_avgabund~master$sb_es+master$lsc_mean, data=master)
sb3.lm = lm(master$sb_avgabund~master$sb_es+master$LAI4_forest+master$lsc_mean, data=master)
sb4.lm = lm(master$sb_avgabund~master$sb_es+master$lc_avgabund, data=master)



#Perform a multiple step wise regression of residual values against HEMI, coastdist, dem80, and slope

cr.lm = lm(master$cr.rq.residuals~master$canopen_road+master$coastdist+master$dem80+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
cr2.lm = lm(master$cr.rq.residuals~master$canopen_road+master$dem80+master$LAI4_road+master$bc05micro+master$bcdiff, data=master)
cr3.lm = lm(master$cr.rq.residuals~master$bc05micro, data=master)

gq.lm = lm(master$gq.rq.residuals~master$MODIS_mean, data=master)
gq2.lm = lm(master$gq.rq.residuals~master$LAI4_road+master$bcdiff+master$logs, data=master)
gq3.lm = lm(master$gq.rq.residuals~master$bc05micro+master$bcdiff, data=master)
gq4.lm = lm(master$gq.rq.residuals~master$bc05micro, data=master)

#lc.lm = lm(master$lc.rq.residuals~master$lsc_mean+master$canopen_road+master$coastdist+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
lc2.lm = lm(master$lc.rq.residuals~master$lsc_mean+master$MODIS_mean+master$bc05micro+master$sb.rq.residuals, data=master)
lc3.lm = lm(master$lc.rq.residuals~master$bc05micro, data=master)

#sb.lm = lm(master$sb.rq.residuals~master$lsc_mean+master$canopen_road+master$coastdist+master$dem80+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
sb2.lm = lm(master$sb.rq.residuals~master$lsc_mean+master$LAI4_forest+master$lc_es, data=master)
#sb3.lm = lm(master$sb.rq.residuals~master$bc05micro, data=master)

sbcomp.lm = lm(master$sb.rq.residuals~master$lc.rq.residuals+master$gq.rq.residuals, data=master)
lccomp.lm = lm(master$lc.rq.residuals~master$sb.rq.residuals+master$gq.rq.residuals, data=master)
gqcomp.lm = lm(master$gq.rq.residuals~master$lc.rq.residuals+master$sb.rq.residuals, data=master)

png(paste(out.dir,"gqrqfittedvses.png", sep=""), height=800, width=800)

plot(master$gq_avgabund, gq.rq$fitted.values,xlab="avgabund", ylab="gq.rq fitted", main="GNYQUEE")
plot(master$gq_avgabund, gq2.rq$fitted.values,xlab="avgabund", ylab="gq.rq2 fitted", main="GNYQUEE")
plot(master$gq_es, gq2.rq$fitted.values,xlab="es", ylab="gq.rq2 fitted", main="GNYQUEE")
plot(master$gq_es, gq.rq$fitted.values,xlab="es", ylab="gq.rq fitted", main="GNYQUEE")

dev.off()

png(paste(out.dir,"sbfitted.png", sep=""), height=800, width=800)
par(mfrow=c(2,2))
plot(master$sb_avgabund, sb.rq$fitted.values,xlab="avgabund", ylab="sb.rq fitted", main="SAPBASI")
plot(master$sb_avgabund, sb2.rq$fitted.values,xlab="avgabund", ylab="sb2.rq fitted", main="SAPBASI")
plot(master$sb_es, sb.rq$fitted.values,xlab="es", ylab="sb.rq fitted", main="SAPBASI")
plot(master$sb_es, sb2.rq$fitted.values,xlab="es", ylab="sb2.rq fitted", main="SAPBASI")

png(paste(out.dir,"physcurves.png", sep=""), height=800, width=800)
par(mfrow=c(2,3))
plot(master$sb_avgabund, sb.rq$fitted.values,xlab="avgabund", ylab="sb.rq fitted", main="SAPBASI")



dev.off()

sb.lm = lm(master$sb_avgabund~master$sb_es, data=master)
sb2.lm = lm(master$sb_avgabund~master$sb_es+master$LAI4_forest+master$lsc_mean+master$MODIS_mean, data=master)
png(paste(out.dir,"sbfittedcompare.png", sep=""), height=800, width=800)
plot(sb2.lm$fitted.values, sb.lm$fitted.values,xlab="OLS2fit", ylab="OLS1fit", main="SAPBASI")
abline(lm.straight, col="red")
dev.off()
png(paste(out.dir,"sblmcompare.png", sep=""), height=500, width=1000)
par(mfrow=c(1,3))
plot(master$sb_es, master$sb_avgabund,xlab="ES", ylab="avgabund", main="SAPBASI")
plot(master$sb_avgabund, sb.lm$fitted.values,xlab="avgabund", ylab="sb.lm fitted", main="SAPBASI")
plot(master$sb_avgabund, sb2.lm$fitted.values,xlab="avgabund", ylab="sb2.lm fitted", main="SAPBASI")
test=matrix(c(0,1,2,3),4,2)
dev.off()


lc.lm = lm(master$lc_avgabund~master$lc_es, data=master)
lc2.lm = lm(master$lc_avgabund~master$lc_es+master$bc05micro+master$lsc_mean, data=master)
png(paste(out.dir,"lcfittedcompare.png", sep=""), height=800, width=800)
plot(lc2.lm$fitted.values, lc.lm$fitted.values,xlab="OLS2fit", ylab="OLS1fit", main="SAPBASI")